/*
 * solver.c
 *
 *  Created on: 30/12/2011
 *      Author: marcio
 */

#include "solver.h"

double*** x;
double* mu;
double* Lx;
int noI;//melhoria
double bestLB;//o melhor limite inferior encontrado
double bestUB;//o melhor limite superior encontrado
double lb;//o limite inferior corrente
double ub;//o limite superior corrente
double pi;//passo do lagrangeano

int compar(const void* a, const void* b){
	int v = *(int*)a;
	int u = *(int*)b;
	return mu[v] - mu[u];
}


double calculateUB(){
	int V[n];
	int cor[n];
	int i,j;
	int color = 1;
	int adj[n];
	int cand = -1;
	for(i = 0; i < n; i++){
		V[i] = i;
		cor[i] = 0;
		adj[n] = 0;
	}
	qsort(V, n, sizeof(int), compar);
	cor[0] = 0;
	for(i = 1; i< n; i++){
		for(j = 0; j < i; j++){
			if(!hasEdge(V[i], V[j])){
				adj[cor[j]] = 1;
			}
		}
		for(j = color - 1; j >= 0; j--){
			if(adj[j] == 0)cand = j;
		}
		if(cand == -1){
			cand = color;
			color++;
		}
		cor[i] = cand;
		cand = -1;
		for(j = 0; j < n; j++){
			adj[j] = 0;
		}
	}
	if(bestUB > k*color){
		bestUB = k*color;
	}
	return k*color;
}

int noImproviment(){
	noI++;
	if(noI > 30){
		pi = pi/2.0;
		noI = 0;
	}
	return 0;
}

double calculateLB(){
	int i;
	double value = 0.0;
	for(i = 0; i < n; i++){
		Lx[va[i]] = solveSingleLP(i);
		value += Lx[va[i]];
		value += k*mu[va[i]];
	}
	if(value > bestLB) bestLB = value;
	else noImproviment();
	return value;
}

int updateMuLambda(){
	double g[nverts()];
	double theta = 0.0;
	double theta_div = 0.0;
	int i,j;
	for(i = 0; i < n; i++){
		g[va[i]] = (double)k;
		for(j = 0; j <= i; j++){
			g[va[i]] -= k*x[va[j]][va[i]];
		}
		theta_div += g[va[i]]*g[va[i]];
	}
	if(theta_div >=  0.001)theta = pi*(ub - lb)/theta_div;
	else theta = pi/n;
	//calcula subgradientes
	for(i = 0; i < n; i++){
		mu[va[i]] += theta*g[va[i]];
		if(mu[va[i]] < 0.0000){
			mu[va[i]] = 0.0;
			if(g[va[i]] <= 0.0) g[va[i]] = 0.0;
		}
	}
	theta_div = 0.0;
	//recalcula subgradente
	for(i = 0; i < n; i++){
		theta_div += g[va[i]]*g[va[i]];
	}
	if(theta_div >=  0.001)theta = pi*((1.05)*ub - lb)/theta_div;
	else theta = pi/n;
	//recalcula subgradientes
	for(i = 0; i < n; i++){
		mu[va[i]] += theta*g[va[i]];
		if(mu[va[i]] < 0.0000){
			mu[va[i]] = 0.0;
		}
	}
	return 0;
}


/*cria a formulação*/
int createFormulation(){
	int i, j;
	pi = 2.0;
	Lx = malloc(sizeof(double)*nverts());
	mu = malloc(sizeof(double)*nverts());
	x = malloc(sizeof(double**)*nverts());
	for(i = 0; i < nverts(); i++){
		mu[i] = 1.0;
		x[i] = malloc(sizeof(double*)*nverts());
		for(j = 0; j < nverts(); j++){
			x[i][j] = malloc(sizeof(double)*k);
		}
	}
	init();
	return 0;
}

/*resolve a formulação*/
double solve(){
	int i;
	while(pi > 0.005){
		changeObjCoef(mu);
		getX(x);
		lb = calculateLB();
		ub = calculateUB();
		updateMuLambda();
	}
	*sol = bestUB;
	return bestLB;
}
